Expresión diferencial

Objetivos

• Conocer el flujo de trabajo empleado en el análisis bioinformático de datos obtenidos por RNA-seq
• Explicar el análisis de expresión diferencial empleando edgeR o DESeq2
• Analizar un conjunto de datos de expresión empleando un script de edgeR

Introducción

El avance en las tecnologías de secuenciación en los últimos años, ha permitido estudiar mayor número de secuencias genéticas provenientes de distintos organismos. Asímismo, los avances en el almacenamiento de la información secuenciada ha favorecido el análisis masivo de datos.

• Genómica (Secuenciación de genomas completos)
• Transcriptómica (Secuenciación de transcriptomas completos)
• Proteómica (Secuenciación de proteomas completos)

El sufijo “omica” refiere el estudio completo de dicho nivel molecular.

Transcriptómica

La secuenciación de transcriptomas completos permite conocer los genes que expresa una célula o un conjunto de células en un momento determinado y bajo condiciones particulares.

Algunas de las preguntas que pueden responderse con el análisis de transcriptomas son:

• Co-expresión de genes entre grupos de muestras y asociación con fenotipos
• Descubrimiento de nuevos genes e isoformas
• Fusión de genes
• Los genes o isoformas diferencialmente expresadas entre dos condiciones (Celúlas A vs Células B, Tratamiento A vs Tratamiento B, Tiempo A vs Tiempo B)

La técnica que permite realizar la secuenciación de los transcriptomas es la secuenciación del RNA (RNA-seq).

Una vez secuenciadas las bibliotecas se genera una cantidad masiva de datos ¿qué a hacer con toda la información obtenida?

Análisis bioinformático de los datos de secuenciación

La cantidad de datos que se generan a partir de la secuenciación de RNA es inmensa. El análisis de expresión diferencial requiere seguir el siguiente pipeline o flujo de trabajo:

pipeline

• 1. Control de calidad de las secuencias crudas
• 2. Filtrado y limpieza de las secuencias En los siguientes vínculos encontrarán los manuales de Trimmomatic y Cutadapt.
• 3. Alineamineto Para una revisión detallada del algoritmo de STAR, pueden revisar el artículo publicado por sus desarrolladores (Dobin et al. 2013). Asimismo la versión más reciente del manual.
• 4. Conteo y estimación de la abundancia. Nuevamente, para una revisión más profunda del algoritmo de RSEM, pueden revisar el siguiente artículo (Li and Dewey 2011) y el manual

Al finalizar el proceso de estimación de la abundancia, RSEM nos devuelve los siguientes archivos por muestra:

## t5.genes.results
## t5.isoforms.results

Y tienen el siguiente contenido de archivos de texto plano:

## gene_id  transcript_id(s)    length  effective_length    expected_count  TPM FPKM
## ENSG00000000003  ENST00000373020,ENST00000494424,ENST00000496771,ENST00000612152,ENST00000614008 2063.22 1892.73 369.00  19.99   14.14
## ENSG00000000005  ENST00000373031,ENST00000485971 940.50  770.06  0.00    0.00    0.00
## ENSG00000000419  ENST00000371582,ENST00000371584,ENST00000371588,ENST00000413082,ENST00000466152,ENST00000494752 1086.32 915.83  853.00  95.51   67.55

Se requieren importar los datos generados (genes.results o .transcripts.results) por RSEM a R. Para ello se requiere de la librería de tximport [(??? ctbTximport2017). Esta librería importa archivos de cuentas generados por:

• Kallisto
• Salmon
• Cufflinks
• Rsem
• … entre otros

Se requiere de una tabla de metadatos:

samples

##    Sample Cell Treatment Replicate     unique_id
## 1      t1    A   Control         1   A_Control_1
## 2      t2    A   Control         2   A_Control_2
## 3      t3    A   Control         3   A_Control_3
## 4      t4    M   Control         1   M_Control_1
## 5      t5    M   Control         2   M_Control_2
## 6      t6    M   Control         3   M_Control_3
## 7      t7    A Verafinib         1 A_Verafinib_1
## 8      t8    A Verafinib         2 A_Verafinib_2
## 9      t9    A Verafinib         3 A_Verafinib_3
## 10    t10    M Verafinib         1 M_Verafinib_1
## 11    t11    M Verafinib         2 M_Verafinib_2
## 12    t12    M Verafinib         3 M_Verafinib_3

Y crear una ruta hacia la ubicación de los archivos:

##Se requiere crear una ruta al directorio en donde se encuentran los archivos
files <- file.path("../RSEM/RSEM/", samples$Sample, paste0(samples$Sample, ".genes.results"))
names(files) <- samples$unique_id
##Se corrobora que la ruta fue creada y contiene los nombres de los archivos correctamente
data.frame(files = files, Exists = file.exists(files))

##                                             files Exists
## A_Control_1     ../RSEM/RSEM//t1/t1.genes.results   TRUE
## A_Control_2     ../RSEM/RSEM//t2/t2.genes.results   TRUE
## A_Control_3     ../RSEM/RSEM//t3/t3.genes.results   TRUE
## M_Control_1     ../RSEM/RSEM//t4/t4.genes.results   TRUE
## M_Control_2     ../RSEM/RSEM//t5/t5.genes.results   TRUE
## M_Control_3     ../RSEM/RSEM//t6/t6.genes.results   TRUE
## A_Verafinib_1   ../RSEM/RSEM//t7/t7.genes.results   TRUE
## A_Verafinib_2   ../RSEM/RSEM//t8/t8.genes.results   TRUE
## A_Verafinib_3   ../RSEM/RSEM//t9/t9.genes.results   TRUE
## M_Verafinib_1 ../RSEM/RSEM//t10/t10.genes.results   TRUE
## M_Verafinib_2 ../RSEM/RSEM//t11/t11.genes.results   TRUE
## M_Verafinib_3 ../RSEM/RSEM//t12/t12.genes.results   TRUE

Cuidado: Revisar bien el orden (jerarquía) y nombre de los directorios y archivos.

Se importan los archivos empleando tximport:

txi.rsem <- tximport(files, type = "rsem", txIn = F, txOut = F)

## reading in files with read_tsv

## 1 2 3 4 5 6 7 8 9 10 11 12

str(txi.rsem)

## List of 4
##  $ abundance          : num [1:60675, 1:12] 16.8 0 118.2 9.8 19.4 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:60675] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457" ...
##   .. ..$ : chr [1:12] "A_Control_1" "A_Control_2" "A_Control_3" "M_Control_1" ...
##  $ counts             : num [1:60675, 1:12] 338 0 1101 319 455 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:60675] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457" ...
##   .. ..$ : chr [1:12] "A_Control_1" "A_Control_2" "A_Control_3" "M_Control_1" ...
##  $ length             : num [1:60675, 1:12] 1981 770 915 3201 2306 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:60675] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457" ...
##   .. ..$ : chr [1:12] "A_Control_1" "A_Control_2" "A_Control_3" "M_Control_1" ...
##  $ countsFromAbundance: chr "no"

##De la lista seleccionamos la matriz de cuentas
counts <- txi.rsem$counts
class(counts)

## [1] "matrix"

head(counts)

##                 A_Control_1 A_Control_2 A_Control_3 M_Control_1 M_Control_2
## ENSG00000000003      338.00      297.00      336.00      492.00      369.00
## ENSG00000000005        0.00        0.00        0.00        0.00        0.00
## ENSG00000000419     1101.00     1199.00     1054.00     1122.00      853.00
## ENSG00000000457      319.26      263.52      342.49      421.71      291.48
## ENSG00000000460      454.74      609.48      416.51      502.29      380.52
## ENSG00000000938        0.00        0.00        0.00        0.00        0.00
##                 M_Control_3 A_Verafinib_1 A_Verafinib_2 A_Verafinib_3
## ENSG00000000003      367.00        327.00        295.00        309.00
## ENSG00000000005        0.00          0.00          0.00          0.00
## ENSG00000000419      917.00        860.00        922.00        889.00
## ENSG00000000457      281.16        361.74        372.45        379.61
## ENSG00000000460      631.84        183.26        190.55        146.39
## ENSG00000000938        1.00          1.00          2.00          1.00
##                 M_Verafinib_1 M_Verafinib_2 M_Verafinib_3
## ENSG00000000003        594.00        609.00         633.0
## ENSG00000000005          1.00          3.00           1.0
## ENSG00000000419        787.00        838.00         762.0
## ENSG00000000457        483.88        492.25         436.6
## ENSG00000000460        216.12        245.75         228.4
## ENSG00000000938          0.00          0.00           0.0

Para el análisis de expresión diferencial requerimos que los valores de las cuentas sean números enteros. Recodificar los valores a integer y convertir la matriz de cuentas a un data frame:

storage.mode(counts) <- "integer"
##Recodificar la matriz de cuentas a data frame
counts <- as.data.frame(counts)
class(counts)

## [1] "data.frame"

names(counts)

##  [1] "A_Control_1"   "A_Control_2"   "A_Control_3"   "M_Control_1"  
##  [5] "M_Control_2"   "M_Control_3"   "A_Verafinib_1" "A_Verafinib_2"
##  [9] "A_Verafinib_3" "M_Verafinib_1" "M_Verafinib_2" "M_Verafinib_3"

Filtrado de cuentas Filtrar los datos es muy útil. Los genes cuyos valores de cuentas son igual a 0 no aportan ninguna información interesante al análisis. Además, eliminar genes con bajos valores de cuentas nos evita obtener resultados falsos positivos.

##Seleccionando genes con al menos 5 cuentas por millón (cpm) en al menos 2 muestras
keep <- rowSums(cpm(counts) >= 5) >=2
table(keep)

## keep
## FALSE  TRUE 
## 47623 13052

##Cortando los datos originales
counts <- counts[keep, ]

Es recomendable no realizar un filtrado astringente, de lo contrario se recuperan solamente genes “housekeeping”.

Tanto para edgeR como DESeq2 requerimos almacenar las cuentas y algunos metadatos en un objeto de forma de lista.

edgeR

Para construir la lista de edgeR requerimos indicar cuáles son los grupos que se compararán y el número de réplicas por grupo:

colnames(counts)

##  [1] "A_Control_1"   "A_Control_2"   "A_Control_3"   "M_Control_1"  
##  [5] "M_Control_2"   "M_Control_3"   "A_Verafinib_1" "A_Verafinib_2"
##  [9] "A_Verafinib_3" "M_Verafinib_1" "M_Verafinib_2" "M_Verafinib_3"

##Usar los nombres de las columnas de la matriz de cuentas para seleccionar los dos grupos experimentales
groups <- factor(sub("..$", "", colnames(counts)))
table(groups)

## groups
##   A_Control A_Verafinib   M_Control M_Verafinib 
##           3           3           3           3

Crear la lista empleando la matriz de cuentas y los grupos:

edgeRlist <- DGEList(counts = counts,
                     group = groups, 
                     genes = rownames(counts))
str(edgeRlist)

## Formal class 'DGEList' [package "edgeR"] with 1 slot
##   ..@ .Data:List of 3
##   .. ..$ : int [1:13052, 1:12] 338 1101 319 454 47 454 542 683 339 3096 ...
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:13052] "ENSG00000000003" "ENSG00000000419" "ENSG00000000457" "ENSG00000000460" ...
##   .. .. .. ..$ : chr [1:12] "A_Control_1" "A_Control_2" "A_Control_3" "M_Control_1" ...
##   .. ..$ :'data.frame':  12 obs. of  3 variables:
##   .. .. ..$ group       : Factor w/ 4 levels "A_Control","A_Verafinib",..: 1 1 1 3 3 3 2 2 2 4 ...
##   .. .. ..$ lib.size    : num [1:12] 13876857 13302895 15560526 18126099 13713631 ...
##   .. .. ..$ norm.factors: num [1:12] 1 1 1 1 1 1 1 1 1 1 ...
##   .. ..$ :'data.frame':  13052 obs. of  1 variable:
##   .. .. ..$ genes: chr [1:13052] "ENSG00000000003" "ENSG00000000419" "ENSG00000000457" "ENSG00000000460" ...

Posteriormente los datos son normalizados empleando el método de los TMM (weighted Trimmed Mean of M-values). Para ello edgeR busca y elimina los valores atípicos (expresión absoluta de muestra como expresión relativa entre muestra) y calcula los factores de normalización. En este artículo encontrarán de manera detallada la explicación de la normalización por TMM (Evans, Hardin, and Stoebel 2018).

edgeRlist <- calcNormFactors(edgeRlist, method = "TMM")
edgeRlist$samples

##                     group lib.size norm.factors
## A_Control_1     A_Control 13876857    0.9752559
## A_Control_2     A_Control 13302895    0.9746757
## A_Control_3     A_Control 15560526    0.9713069
## M_Control_1     M_Control 18126099    0.9686802
## M_Control_2     M_Control 13713631    0.9601787
## M_Control_3     M_Control 15051090    0.9892419
## A_Verafinib_1 A_Verafinib 13422693    1.0372465
## A_Verafinib_2 A_Verafinib 13914127    1.0222470
## A_Verafinib_3 A_Verafinib 13174976    1.0294891
## M_Verafinib_1 M_Verafinib 14548005    1.0436064
## M_Verafinib_2 M_Verafinib 14587258    1.0356998
## M_Verafinib_3 M_Verafinib 14879963    0.9976993

Es importante inspeccionar los datos y verificar que estén correctamente normalizados. Para ello, se grafica la expresión absoluta vs expresión relativa para cada muestra:

La consistencia de las réplicas la podemos verificar mediante un análisis de componentes principales (PCA) o calculando la correlación (Pearson) que existe entre las muestras:

p <- cpm(edgeRlist$counts, log = T)
p <- pca(p)

biplot(p, lab = colnames(edgeRlist$counts), pointSize = 10,
       title = "PCA")

cormat <- cor(cpm(edgeRlist$counts, log = T))
pheatmap(cormat, border_color = NA, main = "Correlation of replicates")

Sin embargo, cuando la variación entre las réplicas es muy grande debido a “batch effects” es recomendable emplear algoritmos que permitan corregir estas variaciones sin causar grandes modificaciones a las cuentas y que dichos algoritmos consideren la distribución de los datos. Recientemente, se ha descrito una herramienta que corrige esta variación tomando en cuenta la distribución bi-nomial negativa de las cuentas (Zhang, Parmigiani, and Johnson 2020).

Para el análisis de expresión diferencial requerimos:

• Crear una matriz del diseño experimental
• Calcular la dispersión de los datos
• Crear una matriz de contrastes

Matriz de diseño experimental

En esta matriz le vamos a indicar a edgeR cuáles muestras corresponden a un grupo o condición experimental. Dado que algunos experimentos pueden tener diseños muy complejos (una muestra pertenece a un tipo celular y a un tratamiento) empleamos la función interna de R, modelado de matrices:

design <- model.matrix(~0+edgeRlist$samples$group)
##El término ~0 le indica a la función no incluir una columna de intersecciones y solamente incluir tantas columnas como grupos en nuestro diseño experimental
colnames(design) <- levels(edgeRlist$samples$group)

design

##    A_Control A_Verafinib M_Control M_Verafinib
## 1          1           0         0           0
## 2          1           0         0           0
## 3          1           0         0           0
## 4          0           0         1           0
## 5          0           0         1           0
## 6          0           0         1           0
## 7          0           1         0           0
## 8          0           1         0           0
## 9          0           1         0           0
## 10         0           0         0           1
## 11         0           0         0           1
## 12         0           0         0           1
## attr(,"assign")
## [1] 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$`edgeRlist$samples$group`
## [1] "contr.treatment"

Dispersión de los datos

Como edgeR ajusta los datos a un modelo de distribución bi-nomial negativo (parecido a Poisson) se requiere calcular un parámetro adicional de dispersión de los datos.

edgeRlist <- estimateDisp(edgeRlist, design = design, robust = T)
plotBCV(edgeRlist)

edgeR calcula la dispersión a tres niveles:

• Common. Valor representativo de todos los genes
• Trended. Calculada por rangos de nivel de expresión (expresión baja - expresión alta)
• Tagwise. Valor individual para cada gen

Matriz de contrastes

edgeR en los análisis de expresión diferencial puede obviar las comparaciones entre los grupos experimentales. Sin embargo, es bueno contar con una matriz de contrastes para indicarle a edgeR cuáles son las comparaciones que queremos hacer entre los datos.

contrast <- makeContrasts(
  "Cells" = "A_Verafinib - A_Control", 
  levels = edgeRlist$design
)

contrast

##              Contrasts
## Levels        Cells
##   A_Control      -1
##   A_Verafinib     1
##   M_Control       0
##   M_Verafinib     0

En este caso, solo tenemos dos condiciones experimentales “Células A” y “Células B”. Con esta matriz de contrastes indicamos que la comparación será buscando los genes diferencialmente expresados en B con respecto a A. Sin embargo, si tuvieramos mayor cantidad de grupos experimentales los podemos comparar, por ejemplo:

Células CT = X, Y, Z

Células Tx = P, Q, R

makeContrast("Cells" = "(P + Q + R)/3 - (X, Y, Z)/3")

El objetivo es que los coeficientes sumen 0

Análisis de expresión diferencial

Los datos deben ajustarse a un modelo lineal bi-nomial negativo. Para ello se utilizará la función glmQLfit con la cual se tiene un control más robusto del “error rate”.

fit <- glmQLFit(edgeRlist, design = design, robust = T, dispersion = edgeRlist$trended.dispersion)

Se realiza la comparación para obtener los genes diferencialmente expresados entre una condición y otra. Con la función glmQLFTest, estamos probando la hipótesis nula “el valor de |lfc| del gen A es igual a 0”. Por lo tanto, los valores de p y p.adj calculados están hechos con respecto a dicha hipótesis nula.

qlf.BvsA <- glmQLFTest(fit, ##Objeto en forma de lista con los datos ajustados a un modelo bi-nomial negativo
                       contrast = contrast[, "Cells"]) ##Matriz de contrastes
##Obtenemos los DEG con expresión distinta de 1, valor de p menor a 0.05, corrigiendo el valor de p por el método de Benjamini-Hochberg
deg.BvsA <- decideTestsDGE(qlf.BvsA, p.value = 0.05, adjust.method = "BH", lfc = 0)
table(deg.BvsA)

## deg.BvsA
##   -1    0    1 
## 3584 5918 3550

##Guardamos los resultados en un data.frame
DEG.BvsA <- DEGResults(qlf.BvsA)
##Añadimos una columna extra a los resultados anteriores indicando la condición de expresión respecto al lfc y al FDR
DEG.BvsA <- edgeResults(DEG.BvsA, logfc = 0, padj = 0.05)

Los datos los visualizamos graficando un “Volcano plot”

volcano_edgeR(DEG.BvsA, lfc = 0, padj = 0.05)

Un procedimiento muy común es seleccionar aquellos genes cuya expresión difeerencial haya sido significativa y que cumplan con un valor de lfc. Por ejemplo: |lfc| > 1 (es decir genes con cuya expresión es > 2 o < 0.5 respecto al CT)

volcano_edgeR(DEG.BvsA, lfc = 1, padj = 0.05)

significant.1 <- DEG.BvsA %>% filter(logFC > 1 & FDR < 0.05 | logFC < -1 & FDR < 0.05)
paste(length(significant.1$genes), "is the number of significant genes")

## [1] "3161 is the number of significant genes"

¡Cuidado! Filtrar los resultados de esta manera es incorrecto. Recordar que en el análisis de expresión diferencial probamos la hipótesis nula “El |lfc| del gen A es igual 0” y los valores de p y FDR corresponden a dicha prueba. El asumir que estos genes filtrados tienen |lfc| > 1 con el valor de FDR previamente calculado provoca un sesgo. Favorece genes con baja expresión y alta variabilidad, invalidando el poder estadístico de la prueba. En otras palabras, incrementamos el riesgo de captar resultados falsos positivos. Aquí encontrarán un artículo en donde se explica detalladamente las implicaciones de realizar este tipo erroneo de filtrados (Love, Huber, and Anders 2014).

Para resolver este problema, tanto edgeR como DESeq2 implementan funciones en donde se prueba la hipótesis nula “el |lfc| del gen A es distinto a x”:

qlf.BvsA.lfc1 <- glmTreat(fit, 
              contrast = contrast[, "Cells"], 
              lfc = 1)
deg.BvsA.lfc1 <- decideTestsDGE(qlf.BvsA.lfc1, p.value = 0.05, adjust.method = "BH", lfc = 1)
table(deg.BvsA.lfc1)

## deg.BvsA.lfc1
##    -1     0     1 
##   792 11343   917

De esta manera podemos seleccionar los genes que estadísticamente tienen |lfc| > 1 y robustecer nuestros resultados.

DEG.BvsA.lfc1 <- DEGResults(qlf.BvsA.lfc1)
DEG.BvsA.lfc1 <- edgeResults(DEG.BvsA.lfc1, logfc = 1, padj = 0.05)

Visualizamos estos datos en un nuevo volcano plot

volcano_edgeR(DEG.BvsA.lfc1, lfc = 1, padj = 0.05)

significant.genes <- DEG.BvsA.lfc1 %>% filter(FDR < 0.05 & logFC > 1 | FDR < 0.05 & logFC < -1)
paste("The number of significant genes with |lfc| > 1 is", length(significant.genes$genes))

## [1] "The number of significant genes with |lfc| > 1 is 1709"

Finalmente, los los genes diferencialmente expresados podemos emplearlos para crear mapas de calor “heatmaps”:

##Obtener los nombres o ids de los genes con expresión diferencial
significant.ids <- significant.genes$genes
##Crear una matriz de cuentas normalizadas por cpm empleando las cuentas que se encuentran guardadas en el objeto edgeRlist
significant.cpm <- cpm(edgeRlist$counts, log = T)
##Cortamos los genes con expresión significativa
significant.cpm <- significant.cpm[significant.ids, c(1, 2, 3, 7, 8, 9)]

##Generamos una paleta de colores
RedBlackGreen <- maPalette(low = "green", high = "red", mid = "black")
##Obtenemos el heatmap
pheatmap(significant.cpm, 
         border_color = NA, 
         color = RedBlackGreen, 
         show_rownames = F, 
         scale = "row", 
         angle_col = 0)

Una vez obtenida la lista de genes con expresión diferencial significativa se procede a realizar análisis de representación de vías (ORA o GSEA) para darle una explicación biloógica a los resultados obtenidos.

DESeq2

En DESeq2 requerímos crear una tabla de metadatos (distinta a la que empleada para importar los datos) con la siguiente información:

##Crear vectores indicando el nombre de las muestras (líneas celulares) y el número de réplicas en cada condición
cell_lines <- c(rep("A", 6))
##Crear un vector indicando las dos condiciones experimentales a comparar
condition <- c(rep("control", 3), rep("tx", 3))

##En un df alojar los datos de las líneas celulares y las condiciones experimentales a las que pertenecen
meta_cells <- data.frame(cell_lines, condition)
##Especificar el nombre de las filas de acuerdo al nombre de las columnas de la matriz de cuentas
rownames(meta_cells) <- names(counts_DESeq)
meta_cells

##               cell_lines condition
## A_Control_1            A   control
## A_Control_2            A   control
## A_Control_3            A   control
## A_Verafinib_1          A        tx
## A_Verafinib_2          A        tx
## A_Verafinib_3          A        tx

Al igual que en edgeR, DESeq2 guarda o aloja los datos de las cuentas y metadatos en un objeto con clase de lista:

DDS_list <- DESeqDataSetFromMatrix(countData = counts_DESeq, ## Datos de las cuentas crudas
                                   colData = meta_cells, ## Tabla de metadatos
                                   design = ~condition) ## Grupos experimentales a los que pertence cada muestra

En este punto tenemos que especificarle a DESeq cuál es el grupo de referencia para realizar las comparaciones:

DDS_list$condition <- relevel(DDS_list$condition, ref = "control")

Para filtrar los genes con expresión baja, DESeq emplea una función interna de normalización para que los datos de las cuentas sean comparables entre ellos:

##Usamos los mismos criterios que en edgeR
filtered <- rowSums(counts(DDS_list) >= 5) >= 2
table(filtered)

## filtered
## FALSE  TRUE 
## 41753 18922

DDS_list <- DDS_list[filtered, ]

Las muestras son normalizadas empleando la función DESeq. Esta función realiza tres pasos en un solo comando:

• Estimación de los size factors (Parecido a los TMMs)
• Estimación de la dispersión
• Ajuste de los datos a un modelo bi-nomial negativo

DDS_list <- DESeq(DDS_list)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

Finalmente, obtenemos los genes diferencialmente expresados. En este caso vamos a probar la hipótesis nula “El lfc del gen A es igual a 0”:

DEG.VerafinibvsCT <- results(DDS_list, 
                         alpha = 0.05, 
                         lfcThreshold = 0)

summary(DEG.VerafinibvsCT)

## 
## out of 18922 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 4029, 21%
## LFC < 0 (down)     : 4162, 22%
## outliers [1]       : 8, 0.042%
## low counts [2]     : 0, 0%
## (mean count < 2)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

Visualizamos los datos en un volcano plot:

##convertimos los resultados a un df
DEG.VerafinibvsCT <- as.data.frame(DEG.VerafinibvsCT) %>%
  rownames_to_column(var = "ensgene_id")
DEG.VerafinibvsCT <- DESEqResults(DEG.VerafinibvsCT, logfc = 0, FDR = 0.05)
volcano_DESeq2(DEG.VerafinibvsCT, lfc = 0, FDR = 0.05)

Y de manera similar a edgeR, generamos un heatmap con los genes diferencialmente expresados:

significant.genes.2 <- DEG.VerafinibvsCT %>% filter(log2FoldChange > 0 & padj < 0.05 | log2FoldChange < 0 & padj < 0.05)
significant.ids.2 <- significant.genes.2$ensgene_id
significant.counts <- counts(DDS_list, normalized = T)
significant.counts <- significant.counts[significant.ids.2, ]
pheatmap(significant.counts, 
         border_color = NA, 
         color = RedBlackGreen, 
         show_rownames = F, 
         scale = "row",
         angle_col = 0)

Análisis posteriores

La lista de genes diferencialmente expresados obtenida por edgeR o DESeq2 no aportan información sobre los procesos biológicos que pueden verse afectados en las células tratadas comparado con el control. Para ello, se requiere realizar análisis que le den un contexto biológico a las listas obtenidas. Por ejemplo los análisis de vías son muy útiles para contextualizar los resultados:

• Over-representation analysis (ORA)
• Gene set enrichment analysis (GSEA)
• Protein-protein interactions (PPI)

En este artículo encontrarán la información detallada de los análisis de vías, sus fortalezas, debilidades y alcances de cada método (García-Campos, Espinal-Enríquez, and Hernández-Lemus 2015).

Referencias

Dobin, Alexander, Carrie A. Davis, Felix Schlesinger, Jorg Drenkow, Chris Zaleski, Sonali Jha, Philippe Batut, Mark Chaisson, and Thomas R. Gingeras. 2013. “STAR: Ultrafast Universal RNA-Seq Aligner.” Bioinformatics 29 (1): 15–21. https://doi.org/10.1093/bioinformatics/bts635.

Evans, Ciaran, Johanna Hardin, and Daniel M. Stoebel. 2018. “Selecting Between-Sample RNA-Seq Normalization Methods from the Perspective of Their Assumptions.” Briefings in Bioinformatics 19 (5): 776–92. https://doi.org/10.1093/bib/bbx008.

García-Campos, Miguel A., Jesús Espinal-Enríquez, and Enrique Hernández-Lemus. 2015. “Pathway Analysis: State of the Art.” Frontiers in Physiology 6: 383. https://doi.org/10.3389/fphys.2015.00383.

Li, Bo, and Colin N. Dewey. 2011. “RSEM: Accurate Transcript Quantification from RNA-Seq Data with or Without a Reference Genome.” BMC Bioinformatics 12 (August): 323. https://doi.org/10.1186/1471-2105-12-323.

Love, Michael I, Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2.” Genome Biology 15 (12): 550. https://doi.org/10.1186/s13059-014-0550-8.

Zhang, Yuqing, Giovanni Parmigiani, and W Evan Johnson. 2020. “ComBat-Seq: Batch Effect Adjustment for RNA-Seq Count Data.” NAR Genomics and Bioinformatics 2 (3): lqaa078. https://doi.org/10.1093/nargab/lqaa078.